# title: greene2005 true random effects estimation algorithm
# author: StableGenius
# date: 2020-01-22
# description: 使用maxLik复现Greene2005提出的真实随机效应模型估计算法
# references:
## - GREENE W. Fixed and Random Effects in Stochastic Frontier Models[J]. Journal of Productivity Analysis, 2005, 23(1): 7–32.
## - GREENE W. Reconsidering heterogeneity in panel data estimators of the stochastic frontier model[J]. Journal of Econometrics, 2005, 126(2): 269–303.

# import packages
library(plm)
library(tidyverse)
library(maxLik)

# Halton draw counts
R <- 200

# calculate log likelihood
cal_logL <- function(epsilon_matrix, dimN, dimT){
  if (!"matrix" %in% class(epsilon_matrix)){
    stop("need matrix!")
  } else if (nrow(epsilon_matrix) != dimN * dimT){
    stop("wrong dimension!")
  }
  
  result <- 0
  for (index in 1:dimN){
    tmp <- epsilon_matrix[((index-1)*dimT + 1) : (dimT*index), ]
    single_result <- apply(tmp, 2, prod) %>%
      sum %>%
      log
    result <- result + single_result
  }
  
  result
}

# "true" random-effect model estimation
sfa_tre <- function(formula, data, model = "product"){
  
  # check model specification: production or cost
  if (model == "product" | model == "production"){
    S <- 1
  } else if (model == "cost"){
    S <- -1
  } else {
    stop("Wrong model specification!")
  }
  
  # check panel data
  if (!"pdata.frame" %in% class(data)){
    stop("data set is not panel data!")
  } else if (!pdim(data)$balanced) {
    stop("unbalanced panel data!")
  }
  
  # dimension of panel data
  # balanced panel only
  dimN <- pdim(data)$nT$n
  dimT <- pdim(data)$nT$T
  
  ## halton draw
  halton_sample <- qnorm(randtoolbox::halton(dimN*dimT, dim = R))
  
  # get variables
  y <- model.response(model.frame(formula, data))
  X <- model.matrix(formula, data)
  
  # parameters initialization
  re_model <- plm(formula = formula, data = data, effect = "individual", model = "random")

  ## beta(with intercept)
  beta_0 <- coef(re_model)
  K <- length(beta_0)
  
  ## sigma
  sigma_0 <- sd(residuals(re_model))
  
  ## lambda
  lambda_0 <- 1
  
  ## sigma_w
  sigma_w <- 1
  
  # parameter initial value (construct and rename)
  parameters <- c(beta_0, sigma_0, lambda_0, sigma_w) 
  names(parameters) <- c(names(beta_0), "sigma", "lambda", "sigma_w")
  
  # log likelihood 
  logL <- function(param){
    param_beta <- param[1:K]
    param_sigma <- param[K+1]
    param_lambda <- param[K+2]
    param_w <- param[K+3]
    
    # restrictive condition
    if(param_sigma<0 | param_lambda<0 | param_w<0)
      return(NA)
    
    # generate random effects
    w <- param_w * halton_sample
    
    # epsilon matrix
    ep <- as.vector(y - X %*% param_beta) - w
    
    # log likelihood value
    singleL <- 2/param_sigma * dnorm(ep/param_sigma) * pnorm(-S*param_lambda*ep/param_sigma)
    cal_logL(singleL, dimN = dimN, dimT = dimT) - dimN * log(R)
  }

  logL_grad <- function(){
    # do nothing
  }
  
  maxLik(logLik = logL, 
         #grad = logL_grad, 
         start = parameters, 
         method = "BFGS")
}